function [S_i,t,COE]=DetermineResonanceOrbit(DU,N,M,J_target,GM_central)
options=odeset('RelTol',1e-6,'AbsTol',1e-9);

am = DU;
%[] Semi-major axis of the planetary moon

T1 = N*2*pi*sqrt(am^3/GM_central);

as =((N/M*sqrt(am^3))^(2))^(1/3);
% T2 = M*2*pi*sqrt(as^3/GM_central)

a_til = as/am;
%[] normalized semi-major aixs

e = sqrt(1-1/(4*a_til)*(J_target-1/a_til)^2);

COE = [as;e;0;0;0;0];

So = Coe2State(COE,GM_central);

[t,S_i] = ode45(@(t,S)R2bpCartesianModel(t,S,GM_central),[0,N*T1],So,options);
S_i = S_i';

end

